Probabilistic ballistic annihilation with continuous velocity distributions 
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We investigate the problem of ballistically controlled reactions where particles either annihilate 
upon collision with probability p, or undergo an elastic shock with probability 1 — p. Restricting 
to homogeneous systems, we provide in the scaling regime that emerges in the long time limit, 
analytical expressions for the exponents describing the time decay of the density and the root- 
mean-square velocity, as continuous functions of the probability p and of a parameter related to the 
dissipation of energy. We work at the level of molecular chaos (non-linear Boltzmann equation), 
and using a systematic Sonine polynomials expansion of the velocity distribution, we obtain in 
arbitrary dimension the first non-Gaussian correction and the corresponding expressions for the 
decay exponents. We implement Monte-Carlo simulations in two dimensions, that are in excellent 
agreement with our analytical predictions. For p < 1, numerical simulations lead to conjecture 
that unlike for pure annihilation (p — 1), the velocity distribution becomes universal, i.e. does not 
depend on the initial conditions. 

I. INTRODUCTION 

We consider an assembly of particles that move freely in d-dimensional space between collisions, where only two 
body collisions are taken into account. The purpose of this paper is to present a model that unifies both the dynamics 
of annihilation 

H H II and of hard-sphere gases m using a continuous parameter p e [0,1], the probability 



that two particles annihilate when they touch each other |2J. In the limiting case p = 1, we recover pure annihilation 
dynamics, and for p — the system of hard spheres. In our system in the limit p — > 0, p > (denoted p — > + ), a 
particle will collide elastically many times before being annihilated. Thus the particles have a diffusing-like motion 
before annihilating. 

Another extensively studied class of problems is the one of diffusion-limited annihilation in which diffusing particles 
annihilate on contact with a given rate H, U The simplest case corresponds to the reaction A + A — > 0. The 
number of particles decays, in the long time regime, as a power law n(t) ~ i~S The decay exponent can be exactly 
computed [llj and is £ = min(l,c?/2), where d is the dimension of the system. However, the time decay exponents 
for the density found in our case when p — > + are different from the exponents found in diffusion-limited systems. 
The reason for this difference is that the underlying microscopic mechanisms responsible for diffusion are different. In 
our case, particles which have a bigger velocity modulus have a bigger annihilation rate than the slow particles. The 
velocity dependence of the annihilation rate is not present in the usual diffusion-limited annihilation. 

It was recently shown |5j that in the long time limit, the annihilation dynamics for dimensions higher than one is 
adequately described by the nonlinear Boltzmann equation. This may be understood in a qualitative way by the fact 
that the density of the gas decays as a function of time, so that the packing fraction (which is the total volume occupied 
by the particles divided by the total volume of the system) decreases, and the role played by correlations (re-collisions) 
becomes neglectible. The Boltzmann equation thus becomes relevant at late times. With this phenomenology in mind, 
we conjecture that in the case of probabilistic ballistic annihilation, the Boltzmann equation adequately describes the 
dynamics for p > 0. For p — 0, the resultin g el astic hard sphere system would be correctly described by Boltzmann's 
equation in the low density regime only |r3.ll2|. 

The paper is organized as follows: in Sec. ^ we first introduce the Boltzmann kinetic equation describing the 
probabilistic annihilation dynamics of a homogeneous system in the scaling regime, that corresponds to asymptotically 
large times. We then provide analytical expressions for the exponents £ and 7 governing the algebraic time decay of the 
particle density and the root mean-square-velocity respectively. Next, we give the first non-Gaussian correction 02 to 
the rescaled velocity distribution by means of a Sonine polynomial expansion. This allows to give explicit expressions 
for the exponents £ and 7 up to the first correction in a^. Sect. IIIII shows the results of direct Monte-Carlo simulations 
(DSMC) that are in very good agreement with the analytical results. In the insight of those simulations we clarify the 
ambiguities following from the analytical computation of a-i |l3j and select the simplest and most accurate relation for 
a2- It is numerically shown that unlike for pure annihilation, the first Sonine correction for < p < 1 does not depend 
on the parameter /i characterizing the initial distribution / for small velocities: \im\ v \_> f(v;t — 0) cx |v| M . We 
also show analytical and numerical evidence that the conjecture put forward in p| according to which the exponent 
£ = Ad/ {Ad + 1) becomes exact in the limiting case p — ► + [Jj. Finally, Sect. IIVI contains our conclusions. 
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II. BOLTZMANN KINETIC EQUATION 
A. Scaling regime 

We consider a system made of spheres of diameter a moving ballistically in <i-dimensional space. If two particles 
touch each other, they annihilate with probability p and thus disappear from the system. With probability 1 — p, 
they undergo an elastic collision. The precollisional velocities v** and the postcollisional ones Vj are related in the 
latter case by 

v** = Vi - (V12 • S-)er, (la) 
v 2 * = v 2 + (vi2 • a-)a, (lb) 

where V12 = vi — v 2 is the relative velocity of two particles, and a a unit vector joining the centers of the grains. 
We consider only two body collisions. The initial spatial distribution of particles is supposed to be and assumed to 
remain homogeneous. 

Let /i(vi;f) be the instantaneous single particle distribution function in R d . The Boltzmann equation for our 
homogeneous system free of forcing reads |5J 



^/i(vi! t) = -po d - x J dw 1 dv 2 d^9{a ■ v 12 ) (a ■ v 12 ) /i(v l5 t)/i(v 2 ; i) + (1 - p) 7 c (/i, fx) 



(2) 



where 9 is the Heaviside distribution, v 12 = v 12 /u 12 , and V\ 2 = |v 12 |. The integration with respect to da runs over 
the solid angle. The first term on the right-hand-side of Eq. J2J) describes the annihilation dynamics, and the second 
one the elastic shocks: the collision term I c reads 

Io(hJx) = a^ X f dw 1 dw 2 da6{a-w 12 ){a-w 12 ) [h<yt*\t)fi(y?;t) - /i(vi;i)/i(v a ;t)]. (3) 

We are searching for an isotropic scaling solution of the homogeneous system, where the time dependence of the 
distribution function is absorbed into the particles density n(t) and in the typical velocity v(t) = y/2(v 2 )/d, where 
(v 2 ) is the mean squared velocity. This imposes the scaling form [HI Il4j 



v d (t)- 



/i(v;t) = =^/(c), (4) 



where the rescaled velocity is given by c = v/v(t). By construction, J f = 1. 



dn 

— = -pu)(t)n, 



B. Decay exponents in the scaling regime 

By integrating Eq. J2J over vi with weights 1 and v 2 , we obtain the number density and energy time evolution 

(5a) 
(5b) 

cm. 

where the collision frequency to is given by 

(6) 



d(nv 
dt 



-paoj(t)nv 2 



uj(t) = n(t)v(t)a d 1 / dc 1 dc 2 da- {a ■ c 12 ) 9 (<r ■ c 12 ) /(ci)/(c 2 ) 

and the energy dissipation parameter a by 

_ J dc x dc 2 da (<r • c 12 ) 9 (S ■ c 12 ) cj/(c 1 )/(c 2 ) 



/ dcc 2 f(c) J dcxdc 2 da- (a ■ c 12 ) 9 (a ■ c 12 ) /(ci)/(c 2 ) 



(7) 



We made use of the fact that the elastic dynamics does not contribute to the decay of energy or density, thus the 
integration over the elastic collision term vanishes. The resolution of Eqs. © follows the method of Ref. [5| and we 
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obtain 



n 
n 

v 

VQ 



l+p- 



1 + a 



-wot 



l+p- 



1 + a 



-U) t 



-2/(l+a) 



(l-a)/(l+a) 



(8a) 
(8b) 



where luq = u>(t = 0) and vq = v(t = 0). We conclude from this result that the dynamics are up to the time rescaling 
t — > t/p (and importantly up to the numerical value of a) the same as the ones obtained for pure annihilation 
The decay exponents are given by n{t) oc and v(t) oc t r , with 
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1 + a' 
a- 1 
a+1' 



The scaling exponents consequently satisfy the constraint £ + 7=1. 



(9a) 
(9b) 



C. Rescaled kinetic equation 

Inserting the scaling form Q in Eq. J2J and making use of Eqs. 1(8} , we obtain after some algebra 



C12 



./(>, ) ^ ^ /(.n i / dc 2 | Cl2 |/(c 2 ) - 



where 



and 



HfJ)= dc 2 / dS0(S-c 12 )(S-c 12 ) /(cD/(cr)-/(ci)/(c 2 ) 



ft 



/ da 6»(?-vi 2 )(S ; -9i 2 ) = -— — - 

j R d r[(d + i 



)/2]' 



(10) 



(11) 



(12) 



r being the gamma function. In equation (|10|) . the angular brackets denote average with weight /: for a given function 
<?(ci,c 2 ) 



(q) = / dcidc 2 g(ci,c 2 )/(ci,c 2 ) 



Making use of the identity [T^ 



dcc fc (d + c^- ) /(c) = -k(c k ), 



and integrating Eq. 11U|) over C! with weight e[ , one obtains 

a 2 f (c 12C fe) ^ , 1-p 2 



Vfc > 0, 



where = — J" Rd dcicf /(/,/) and a = (ci2cf)/((ci2)(cf)) is the energy dissipation parameter. 



(13) 



(14) 



(15) 



D. First non-Gaussian correction 



The solution of the Boltzmann equation for pure annihilation dynamics (p = 1) is non Gaussian in several aspects. 
The tail of the distribution is overpopulated |5j , and deviations from the Gaussian behavior may also be observed 
near to the velocity origin 0, @ ■ It is thus reasonable to think that the velocity distribution function obtained upon 
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solving Eq. (|10fl will show similar deviations. To study the correction close to the origin, it is convenient to apply a 
Sonine expansion for the distribution function /(c) |Tfij 



m = M(c) 



i>l 



aiSi(c 2 ) 



(16) 



where M.{c) — n d l 2 exp(— c 2 ) is the Maxwellian, and Si(c 2 ) the Sonine polynomials. Due to the constraint 



= d/2, 

the first correction a\ vanishes |14| . For our purposes, it is sufficient to push the truncation of expression i|lt)|) to second 
order, where S 2 {x) = x 2 /2 — (d + 2)x/2 + d(d + 2) /8. In order to compute a and a 2 , one may follow the method used 
for inelastic granular gases in Ref. UJ : we may use the hierarchy (|15[) for k = 2 and k — 4 to obtain a system of two 
equations for the two unknowns a and a 2 . The calculations are however tedious and it appears useful to consider the 
alternative method that consists in invoking the limit of vanishing velocities of Eq. i|10|) 5] . Indeed, since we expect 
that the tail of the exact solution for the distribution function differs significantly from .M(c)[l + J2i>i a iSi{c 2 )\, the 

computation of low order moments of / should give a more accurate result. From Eq. Q 



1 



,1 — a 



/(0) = /(0)( Cl ) - 



1-pl 



(17) 



We see that the limit in Eq. H17|) involves moments of a lower order than 114. Neglecting the corrections a», i > 3, the 
computation of the latter limit gives (see Appendix) 



hm /(/,/) = 
ci-o U,JJ 20F 



-a 2 



d(d 



l(i 



(18) 



where Sd — 2it d / 2 /T{d/2) is the surface of the d-dimensional sphere. Inserting Eq. (|18fl in Eq. (|17|) . one obtains a 
relation between a and a 2 that is supplemented with that corresponding to k = 2 in 115|) . in order to finally obtain 
a and a 2 . To this end, we make use of the various relations between moments of the velocity distribution and the 
fourth cumulant 02 derived in To linear order in a 2 , the corresponding system reads 




(d-1) 



0(a 2 ), 



(19) 
(20) 



where use have been made of the relation /i2 = (the elastic shocks conserve the total kinetic energy of the colli ding 
pairs), which consequently eliminates p in the second relation. However, as it was shown in previous works |l3t Il6| . 
there are some ambiguities arising from the linearization procedure, that may affect a 2 if this quantity is not small 
enough. We have thus solved the full nonlinear problem, and then in order to have a simpler expression of a 2 , chosen 
the linearizing scheme that yields the closest result (the difference does not exceed 10%) to the nonlinear solution. It 
turns out as well that this scheme is the closest one to the numerical simulations of Sect. Hill This correction is given 
by: 



02 (p) = 8 



3- 2\/2 



4d + 6 - V2 + ±j£8V2(d - 1 



(21) 



In the limiting case of pure annihilation p — > 1, one recovers the result of Ref. 

Inserting this result into the definition Eqs. (J5J, we obtain the decay exponents £ and 7 = 1 — £. In the limit p — > + , 
we note that a 2 vanishes, as may have been anticipated: the velocity distribution then becomes close to its elastic 
Maxwellian counterpart that holds for p = 0. In this limit, the decay exponent is £ = 4d/(Ad+ 1), as conjectured in 
0. We emphasize that the limit p — * is singular: £ is bounded from above by 4d/ (4d + 1) for any p > 0, whereas 
£ vanishes for p = 0. It is therefore important to exclude p — from the limit p — > in order to get well behaved 
limiting expressions. 



III. SIMULATION RESULTS 



We implement a direct Monte-Carlo simulation (DSMC) scheme in order to solve the Boltzmann equation. The 
algorithm may briefly be described as follows. We choose at random two different particles {i, j}. If their velocity 
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FIG. 1: First Sonine correction 02 from the analytical estimate [12 U and from DSMC as a function of the annihilation probability, 
for d — 2. The initial number of particles is 5 x 10 6 , and each value is obtained from approximately 10 4 independent runs. The 
results are not sensitive to the initial velocity distribution. However, the convergence process is much faster starting from a 
Gaussian distribution. 



is such that u = v.y • <r > 0, they may collide. Time is subsequently increased by (iV 2 w) _1 , where N is the number 
of remaining particles. With probability p the two particles are then removed from the system, and with probability 
1 — p their velocity is modified according to Eqs. QJ. For more details on the method see 0, 0, 0, [l^. As the 
number of particles decreases, the statistics at late times suffers from enhanced noise. It is thus necessary to average 
over many independent realizations. 

In dimension one, the dynamics of annihilation creates strong correlations between particles [2(- This precludes a 
Boltzmann approach that relies on the molecular chaos assumption. We will thus focus on numerical simulations of 
two-dimensional systems, and we expect the role of correlations to diminish when the dimensionality increases. 



A. First Sonine correction 



Making use of the relation between 02 and the fourth cumulant of the rescaled velocity distribution 0] 

a2= d(^) (c4) - 1 ' (22) 

we show in Fig. ^the numerical values of the first Sonine correction 02 for different values of p. The agreement with 
Eq. I|21|) is good in most cases. 

It turns out that the discrepancy between Eq. H21(l and DSMC is mainly due to the limit method of computing 
a-i . This method yields a very precise distribution / in the relevant region of interest in the framework of a Sonine 
polynomial expansion, namely the small velocity region. On the other hand, it is less accurate in the less interesting 
high velocity region, hence the discrepancy [T3 |. 



B. Decay Exponents 



Plotting the density n/riQ (and the root-mean-squared velocity v/vq) as a function of time t on a log-log plot gives 
the decay exponents (see Fig. [5J. The numerical results are in agreement with the analytical predictions obtained 
from the set of Eqs. i|19[l and l|2U|) that is inserted in Eq. ©. The predicted power-law behavior is observed over 
several decades, as shown by Fig. |3 for p = 0.5. In Fig.0] we show that the scaling relation £ + 7 = 1 is well obeyed 
for all values of p. Such a relation holds in fact independently of the molecular chaos assumption underlying the 
Boltzmann equation. 
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FIG. 2: The decay exponents £ and 7 (inset) in two dimensions, obtained analytically from Eqs. 119H and 12011 that are 
inserted in Eq. @, and from DSMC (symbols). The initial number of particles is 5 x 10 6 , and the number of independent 
runs approximately 100. The values of the exponents are not very sensitive to the probability p. The horizontal line shows the 
Maxwellian analytical prediction to zeroth order in az, i.e. £ and 7 from 1191 1 and @ with 0,2 = 0. 




FIG. 3: The time dependence in two dimensions of n and v (inset) on a logarithmic scale for p = 0.6 and a Gaussian initial 
velocity distribution, showing a clear power law behavior. The straight line is the linear interpolation giving the decay exponent. 
No (N) is the initial (remaining) number of particles. We have denoted vo the initial root-mean-square velocity v. The same 
quantity is denoted v for t > 0. The deviation observed for large times is due to the low number of remaining particles. 



C. Evolution toward the asymptotic distribution 



In order to have a more precise understanding and accuracy check of our results, it is useful to study the velocity 
distribution in the scaling regime. Indeed, the distribution may be adequately described by the Sonine correction 02 at 
late times only. Before the scaling regime is reached, the velocity distribution f(c±) is time-dependent. A very precise 
check consists in studying the evolution of the non-Gaussianities. To this end, it is useful to consider numerically the 
quantity f(ci)/M(ci) = 1 + OaSaCci)- Fig.[S]shows the evolution of f(ci)/M.(ci) for different times corresponding to 
different densities, starting from an initial Gaussian distribution. 

It turns out that both methods of computing ci2 (directly using its definition in terms of the fourth cumulant iL'l'li 
or using f(ci)/M(ci)), are fully compatible numerically. However, the latter method requires much more extensive 
simulations. It is instructive to investigate the evolution toward the asymptotic solution starting from different initial 
distributions, which are characterized by their behavior near the origin. To this extend we define the exponent fi by 
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FIG. 4: Numerical verification of the relation £ + 7 = 1 in two dimensions for different values of p. Note the j/-scale. 
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FIG. 5: Plot of f(a)/M(ci) at different times corresponding to different densities, for p = 0.5. The initial number of particles 
is 2 x 10 7 and there are approximately 10 s independent runs. The initial distribution is Gaussian and thus corresponds to the 
flat curve. The continuous curve is the analytical prediction I + 0S21S2 with 02 given by Eq. 12 H . The inset shows a magnification 
of the small velocities region. 



the behavior /(c) ~ |c| M for c — > 0. Fig. shows the non-Gaussianities of the evolution towards the scaling function 
for an initial distribution characterized by fi = 3, and Fig. for fi — —3/2. 

For both initial distributions fi = 3 and ji — —3/2, the solution is attracted toward a scaling function characterized 
by ji = 0. Hence, there is a qualitative difference between probabilistic annihilation and pure annihilation. Indeed, 
it was shown in a previous work that for pure annihilation /x is conserved yy, and more importantly that \i indexes 
the "universality classes" of this process (two distributions with the same \x are characterized by the same long time 
exponent £). Obviously, adding the effect of elastic collisions in the dynamics of probabilistic annihilation violates 
the conservation of ji. Next, the question is to know whether the asymptotic distribution depends on fi or not. We 



consequently show in Fig. the ratio / (m=0) (ci)// (ai=3) (ci) = (1 + a^" 0) )/(l 



>=3) 



The ratio tends to unity, which implies that a 



(m=0) 



>=3) 



Moreover, we checked that for the negative value 



fi = —3/2, the same conclusion holds. The convergence is however slower due to the divergence of the initial 
distribution near the velocity origin. We thus conjecture that not only the first Sonine coefficient of probabilistic 
annihilation but also the full velocity distribution (and hence, all decay exponents) show an universal property in the 
sense that they do not depend on the initial velocity distribution if < p < 1. This is a nontrivial result since it was 
shown that this is not true in the case of pure annihilation p = 1 Q . 
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FIG. 6: Same as Fig. E|but for an initial distribution such that /j, — 3. 
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FIG. 7: Same as Fig.|S|but for an initial distribution such that fj, = —3/2 and initially 4 x 10 7 particles. 

Finally, in order to clarify the relevance of the scaling function, we studied the fourth cumulant ai as a function of 
Nq/N, for the same parameters as those in Figs. 15171 The result is shown in Fig. EI The fact that ai reaches a plateau 
indicates that the system enters a scaling regime at late times. For pL = —3/2 (Fig.EJ), due to the initial central peak, 
the initial distribution is extremely different from its late time asymptotic counterpart, so that the transient evolution 
takes longer and the plateau regime is only approached. Finally, it may be observed in Fig. [^1 that for the 3 initial 
conditions the fourth cumulants converge to the same value. This is a further illustration of the universal behaviour 
discussed above. 



IV. CONCLUSIONS 

A system made of spherical particles moving freely in d-dimensional space was studied. When two particles collide, 
they annihilate with probability p or undergo an elastic collision with probability 1 — p. We gave empirical arguments 
for the relevance of the Boltzmann description in this system. We obtained analytically the decay exponents of the 
density of particles and of the root-mean-squared velocity in terms of the energy dissipation parameter a. It turns 
out that upon rescaling time according to t — > t/p, p > 0, the formal structure of the decay equations is the same as 
in the case of pure annihilation p — 1. 

In the scaling regime (that emerges in the long time limit), the first Sonine correction ai to the Maxwellian 
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FIG. 8: Plot of / (m " 0) (ci)// i>=3) ( c i) for three different times, and p = 0.5. We see that for late times the ratio of the two 
distributions tends to unity, which leads to conjecture that the first Sonine corrections ai are the same in both cases fi = and 
fi = 3. The results reported here correspond to particularly extensive simulations (note the vertical scale) 
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FIG. 9: Plot of ci2 as a function of the densities No/N for different values of (j,. There are approximately 5 x 10 4 independent 
runs. 

distribution was obtained as a function of the continuous parameter p. This allows to establish an explicit relation 
for the decay exponents. It was shown that in the limit p — > + , the exponent £ governing the decay of particles, 
nit) cx t~^, is given by £ = 4e?/(4d + 1), thereby confirming a conjecture put forward in 

Numerical simulations (DSMC) in two dimensions are in agreement with the analytical correction a<z (p). Moreover, 
the analytical values for the decay exponents obtained from the first correction are in good agreement as well with 
numerics. The relation £ + 7=1 was shown to hold for all values of p. The study of the dynamics of non-Gaussianitics 
embodied in aiSi reveals a qualitative difference with pure annihilation dynamics: the parameter /i describing the 
small velocity behavior of the rescaled distribution is not conserved for probabilistic annihilation when < p < 1. 
Numerical results for different values of fj, leads to conjecture the universality of the rescaled velocity distribution in 
this process (this universality being lost for pure annihilation only, i.e. for p = 1). 
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APPENDIX: CALCULATION OF THE LIMIT ci^OOF THE COLLISION TERM I 



This quantity may be obtained as a particular case of a previous calculation [13J. The result is the following. We 
define the loss term and gain term I g by 

(A. la) 
(A.lb) 

so that lim Cl ^o I(f, f) = h + Ig- Within the framework of the Sonine expansion (|lfc>|> and neglecting the coefficients 
di, i > 3, the calculation of the latter integrals gives: 



k = - lim / dc 2 / da9(a ■ Ci 2 )(ct • Ci 2 )/(ci)/(c 2 ) 
I g = lim/ dc 2 f da0(a-c 12 )(a-c u )f(cT)f(4*), 



Ii = - 



S d M(0) 
2^ 



1 + a 2 



d{d + 2) 



0.1 

1 - -f 



(A.2) 



~ _ S d M(0) 
9 20F 



1 + a 2 



d 2 - 2d + 3 2 3d(d - 2) 
8 + ° 2 64 



(A.3) 



where Sd = 2n d / 2 /T(d/2) is the surface of the d-dimensional sphere, and T the gamma function. Summing Eqs. <|A.2|) 
and <|A. 3|> leads to the result given by Eq. l|18fl . 
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